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Abstract 

The incompatibilities between the initial and boundary data will cause 
singularities at the time-space corners, which in turn adversely affect the 
accuracy of the numerical schemes used to compute the solutions. We 
study the corner singularity issue for nonlinear evolution equations in ID, 
and propose two remedy procedures that effectively recover much of the 
accuracy of the numerical scheme in use. Applications of the remedy 
procedures to the ID viscous Burgers equation, and to the ID nonlinear 
reaction-diffusion equation are presented. The remedy procedures are 
applicable to other nonlinear diffusion equations as well. 

1 Introduction 

It is well known in the mathematics community that smooth boundary and 
initial data do not guarantee smooth solutions to the initial-boundary value 
problems of time-dependent PDEs. Even if the existence and uniqueness of the 
solution arc proved and all the given data are as smooth as desired, then, in 
order for the solution to be smooth near < = 0, it is necessary and sufficient 
that the boundary data and the initial data satisfy an infinite set of so-called 
compatibility conditions. Sec [9, 10, 11, 12, 13, 14, 15, 17], and below. 

When the boundary and initial data fail to satisfy some of the compatibility 
conditions, singularities may occur at the corner of the time and spatial axes. 
In some cases, for simple problems, this is not a critical issue, because the 
singularity is short-lived. 

However in recent years the issue of the compatibility conditions for the ini- 
tial boundary value problems of time-dependent PDEs started to receive atten- 
tion in the numerical simulation community, because larger and more complex 
problems are handled thanks to the ever-growing computing power that is avail- 
able. This produces a need to better understand what happens during the short 
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initial period for certain physical processes. Boyd and Flyer [3] analyzed the 
connection between incompatibility and the rate of convergence of Chebyshev 
spectral series, and discussed the remedy procedure of smoothing the initial con- 
ditions. Flyer and Swarztrauber [8] studied the effect of the incompatibilities 
on the convergence rate of spectral and finite difference methods. Flyer and 
Fornberg [6] proposed a remedy procedure, based on the idea of singular corner 
functions, for the heat equation, and for variable coefficient convection-diffusion 
equations as well. Flyer and Fornberg [7] studied the corner basis functions for 
some dispersive equations. Bicniasz [2] modified the remedy procedure proposed 
by Flyer and Fornberg [6], and applied it to a diffusion- reaction system arising 
from electrochemistry. 

In this article we study, from a numerical point of view, the compatibility 
issue for nonlinear diffusive PDEs. We shall first present our approach in de- 
tails for the classical viscous Burgers equation in dimension one. However our 
approach does not depend on any particular property of the Burgers equation 
other than its diffusiveness. Hence we believe that it applies to other nonlin- 
ear diffusive equations as well. To demonstrate this, we shall also apply the 
approach to a ID nonlinear reaction-diffusion equation. 

For the Burgers equation the correction procedure proposed in [6] cannot 
be expected to work to its full strength, for two reasons. The first is that, for 
a nonlinear PDE, the singular corner functions, which satisfy the equation and 
also display corner singularities, are hard, if not impossible, to find. The second 
reason is that the superposition property does not hold for nonlinear equations, 
which means that even if in some rare cases we find the singular corner functions 
for the nonlinear equation, we cannot separate them from the solution without 
introducing singular terms into the equation. However, we speculate here that 
the singular corner functions, derived from the linearized Burgers equation, can 
be employed to remove the zero order incompatibility (see Section 2.2). What is 
less obvious is that the singular corner functions can also be employed to remove 
higher order singularities (see Section 2.3). However, as has been said, we cannot 
avoid introducing singular terms into the equation. To overcome the difficulty 
associated with singular terms in the equation we choose the Galcrkin finite 
element method (FEM) as the means for constructing the numerical scheme for 
the equation, because the Galerkin FEM is based on the weak formulation of the 
PDE. and hence is potentially more tolerant of the singularities in the equation. 
The numerical results confirm the effectiveness of the correction procedures that 
wc propose. 

Flyer et. al. [6], Bicniasz [2] and the current study all seek the solutions of 
the target equations (a different equation for each study) in the form 

u = v + S, (I.l) 

where 5 is a linear combination of the singular corner functions of the linear 
diffusive equations. The major difference between our approach and the ap- 
proaches of Flyer et. al. in [6] and Bicniasz in [2] lies in the way the correction 
procedure is implemented. Both Flyer et. al. and Bicniasz derived the differen- 
tial equation for the new unknown v and constructed the numerical scheme for 
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this equation directly. We instead combine the correction procedure with the 
appropriate numerical method, the Galerkin FEM in this study, and work with 
the weak formulation of the equation. In addition, Bieniasz considered only the 
zeroth order incompatibility, while our study considers the incompatibilities up 
to the first order. 

The rest of the article is organized as follows. In Section 2 we recall the 
compatibility conditions for the viscous Burgers equation, and describe the cor- 
rection procedures for the incompatibilities between the initial and boundary 
conditions. In Section 3 we present numerical results to verify the effectiveness 
of the proposed correction procedures. In Section 4 we derive the correction 
procedures for the ID nonlinear reaction-diffusion equation. Numerical results 
are also presented. We conclude with Section 5. 

All the equations considered in this article and in the references quoted above 
relate to evolution equations in space dimension one. In an article in preparation 
[4] we will consider higher spatial dimensions which necessitate totally different 
methods. 

2 The numerical scheme and correction proce- 
dures 

2.1 The Viscous Burgers equation 

Wc consider the initial and boundary value problem for the Burgers equation: 

Ut + uUx — vu^,j. = 0, < a; < 1, t > 0, 
< u{Q,t)=gi{t), u{l,t)=g2{t), t > 0, (2.1) 
^ u{x,0) = h{x), < .T < 1, 

where i/ is a small positive parameter representing the viscosity, and gi , 52 and 
h are given real functions, assumed to be as smooth as desired. The exact 
problem (2.1) is known to have a unique solution for all times; this equation 
is indeed similar, but simpler, than the 2 dimensional incompressible Navicr- 
Stokes Equations for which the existence and uniqueness of the solution is known 
(see e.g. [16]). Looking for a weak solution, if h is given in L^{0, 1) and gi, 32 in 
H^O,T), then u exists and is unique in C{[0,T]; L'^(0, 1)) n L'^{0,T; H^{0, 1)). 
If, as we assume here, h, 51, 52 are smooth, then u is smooth up to C°° regularity 
in [0,1] X (0,r]. Now the fact that the interval (0,r] is open at is related 
to the compatibility problem wc are addressing here; even if h,gi,g2 are given 
C°° (in respectively [0,1] and [0,r]). For the (unique) solution to be smooth 
(even C'^, ^C"^) near i = 0, that is in [0, 1] x [0, T], /i, 32 must satisfy certain 
compatibility conditions as described in the references quoted [13, 14, 17]. We 
now make explicit the first and second compatibility conditions for (2.1) which 
guarantee respectively that u is C'',C^, in [0, 1] x [0,r]. 
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The compatibility conditions require that at the corners of the time and 
spatial axes, the derivatives of the solutions computed through the boundary 
conditions be equal to those computed through the Cauchy-Kowalevski rules, 
that is, for the viscous Burgers equation above. 



1^' order: guiO) ^ -h{0)h'{0) + lyh" (0), g2ti0) ^ 'h{l)h' (!) + i^h" (!) 



We note here that the compatibilities at the left and right corners can be 
treated separately, and the method presented below apply to the incompatibil- 
ities at both corners. For this reason, and for the sake of simplicity, we study, 
in this article, the case in which the compatibility conditions at a; = 1 are met, 
at least to a certain order, but those at the left a; = are not. Assuming so, we 



Nonzero ao, ai represent incompatibilities at the left time-space corner. 

2.2 Correction procedure 1 for the zeroth order incom- 
patibility 

We aim to separate the singular part of the solution from the nonsingular part by 
using certain singular corner functions, but this approach will inevitably intro- 
duce singular terms into the equation, because the singular corner functions for 
the nonlinear equation (2.1)^^ is hard, if not impossible, to find, and also because 
the superposition property does not hold for nonlinear equations. We, however, 
speculate that the correction procedure can be employed to remove the zero or- 
der incompatibility, which is the most significant one. To overcome the difficulty 
associated with the singular terms in the equation we choose Galerkin finite el- 
ement method (FEM) as the means for constructing the numerical scheme for 
(2.1), because Galerkin FEM is based on the weak formulation of the PDE, and 
therefore is potentially more tolerant of the singularities in the equation. 

The weak formulation of (2.1)-^ can be formally derived as follows. Let v be 
a continuous function that vanishes at a; = 0, 1. Multiplying (2.1)-^ by v and 
integrating the equation by parts over (0, 1), we obtain 



0*^ order: gi{0) = h{0), 



52(0) = hil) 



let 



ao = ffi(0)-/i(0) 7^0, 

ai = git(O) + h{0)h'{0) - iyh"{0) ^ 0, 



(wt, v) + -(u^, Vx) + v{ux, Vx) = 0, 



(2.2) 



where (/, g) = /p fg Ax for any integrable functions. 
We introduce the corner function: 




(2.3) 
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and we notice that satisfies the heat equation and displays a singularity at 
(0, 0): 

Sot — vSoxx = 0, 
< So(0,t) = l, (2.4) 
5o(x,0) = 0. 

We let N be the number of segments in the interval [0, 1], and h — l/N, 
and let Vh be the finite element space, with basis {(po, <pi, ■ ■ ■ , ^n}- We look 
for a solution of (2.1) in the form 

u{x,t) = aoSoix^t) +Vhix,t), (2.5) 

where Vh{-,t) G Vt for a.c. t> Q. Imposing the boundary conditions (2.1)3, '^"^^ 
noticing (2.4)2, have 

fu/i(0,t) =5i(t)-ao, 
\vh{l,t) ^ g2{t) - aoSQ{\,t). 

Imposing the initial condition (2.1)3, ^^'^ noticing (2.4)3, ^® have 

Vh{x,0) = h{x), Q<x<l. (2.7) 

We observe that the zeroth order incompatibility at the left time-space corner 
has been removed for Vh] indeed 

Vh{Q,0)^9l{Q)~ao = h{Q). (2.8) 

The compatibility conditions at the right corner are not affected because So is 
smooth there. 

We then plug (2.5) into (2.2) and take v — vu ^Vu, with w/i(0) — w/i(l) = 0, 
and we obtain 

[aoSot + Vht, Vh) - ^((ao-So + vuf , vux) + v{aoSox + Vhx, vhx) = 0. (2.9) 
Noticing that 5*0 satisfies the heat equation, we rewrite (2.9) as 

{vht, Vh) - ^(wL ^hx) - {a-oSoVh, Vhx) + ^{vhx, Vhx) = ]^al{Sl, Vhx)- (2.10) 

Coupling (2.10) with the boundary conditions (2.6) and the initial condition 
(2.7) we can find Vh- 

After we find Vh we recover the original solution of (2.1) by (2.5). 

Numerical experiments will be presented in Section 3. For the tests that we 
have done, the errors during the initial period of time is reduced by a magnitude 
of more than one order. 



5 



2.3 Correction procedure 2 for the zeroth and first order 
incompatibilities 

To further improve the accuracy we consider the next (first order) incompati- 
biUty and introduce the second corner function: 



^1 



So{x,t) dr. 



We notice that the function Si satisfies 

Sit — i^Sixx — 0, 
< SiiO,t)=t, 
5i(x,0) =0. 

We look for a solution of (2.1) in the form 

u{x, t) — Q!oS'o(a;, t) + aiSi{x, t) + Vh{x, t), 



(2.11) 



(2.12) 



(2.13) 



where V}i{-,t) € Vh for almost every t > 0. Imposing the boundary conditions 
(2.1)2, ^^'-^ noticing (2.4)2 ^^'^ (2-12)27 we have 



«/i(0,0 = .91 (i) - ao - ait, 

v,i{l,t) = g2{t) - aoSo{l,t) - aiSi{l,t). 



(2.14) 



Imposing the initial condition (2.1)3, ^^^^ noticing (2.4)3 ^'^"^ (2-12)3, have 

Vh{x,0) = h{x), < a; < 1. (2.15) 

As for the zeroth order correction procedure, this correction procedure also 
removes the zeroth order incompatibility for Vh at the left corner; indeed 



Vh{0,Q)^gi{0)-ao^h{0). 



(2.16) 



The singular corner function Si has no effect on the zero order incompatibility, 
but helps to remove the first order incompatibility. To see this, we plug (2.13) 
into the original equation (2.1)^^, and noticing that 5*0 and Si satisfy the heat 
equations (2.4) and (2.5) respectively, we have 



Vht + ( (ao'S'o + aiSi)vh + ^{aoSo + ctiSi)'^ 



+ VhVhx - vvhxx = 0. (2.17) 



We first calculate i;;it(0, 0) using the boundary condition (2.14)-^, 

i;,,t(0,0) =gu(0)-ai. (2.18) 
Then applying Cauchy-Kowalevsky rule to equation (2.17), we have 

i'/.t(0,0) = 

VhVhx+VVhxx ■ (2-19) 



(ao-So + OLiSi)vh + -(aoS'o + ai<5'i)^ 



(0,0) 
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Noticing that S'o(a:, 0) = 5i(x,0) = (see (2.4)3 (2-12)3), '^^^ have 



Vht{0, 0) = -h{0)h40) + uh^40). (2.20) 

The right hand sides of (2.18) and (2.20) are equal by the definition of ai. 

This correction procedure does remove the zeroth and first order incompat- 
ibilities, but it introduces singular terms into the equation. To overcome this 
difficulty we again construct the numerical scheme by Galerkin FEM, for the 
reason already mentioned before. As for the zeroth order correction procedure, 
we plug (2.13) into (2.2) and take v — Vk & Vh, with ■!;;i(0) = Vh{l) = 0, and we 
obtain 

{vht, Vh) - ^(wL Vhx) - {{aoSo + aiSi)vh, Vhx) + 

v{vhx, Vhx) = ^ ((aoS'o + "iS'i)^, Vhx) ■ (2.21) 

We supplement (2.21) with the boundary conditions (2.14) and initial condition 
(2.15), and solve the resulting system for Vh- Then we recover the solution u of 
(2.1) by (2.13). 

Our results showed that, when this procedure is applied, the accuracy of the 
result is further improved. See Section 3.2. 




3 Numerical implementation of the correction 
procedures and the results 

3.1 The numerical schemes 

To unify the presentation of the numerical schemes for the different correction 
procedures we introduce 

if no correction procedure is applied, 
if Correction procedure 1 is applied, (3.1) 
ao'5'o + OiiSi if Correction procedure 2 is applied. 

The boundary conditions for Vh^ (2.6) or (2.14) depending on the correction 
procedure, can be written in one common form: 

{vh{0,t)^g,{t)-S{Q,t), 
\vh{l,t)^g2{t)^S{l,t). 

The equation for tj/i, (2.10) or (2.21), can also be written in one common form, 

{Vht: Vh) ~ ^(wL Vhx) - {SVh, Vhx) + v{vhx, Vhx) = ^(5*^, Vhx), (3.3) 

for every Vh G Vh with Vh{0) = Vh{l) = 0. 
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The incompatibilities between the initial and boundary conditions have a 
more severe effect on higher order schemes than on lower order schemes (see 
[3, 8]). For the basis of the finite element space Vt (introduced in Section 2.2) 
we choose piecewise linear functions, which usually provide second order approx- 
imations to smooth functions. The effectiveness of the correction procedures is 
already evident with the resulting numerical scheme. We leave the endeavor for 
higher order schemes to future work. 

Let N be the number of segments in the interval [0,1], Ax ~ — , and 
Xj ~ jAx for < j < N . Let c/jj be the piecewise linear hat functions: Then 





x„=0 



Figure 3.1: Piecewise linear finite elements 
we write Vh as a linear combination of these basis functions: 



N 



Vhix,t) ^^Vnit)(l)nix), (3.4) 



n=0 

and w„, for < n < iV, are the unknowns. 

Imposing the boundary conditions (3.2) on we obtain 

( vo ^ gi{t) ~ SiO.t), 
\vN = 92it)-Sil,t). 

Next we plug (3.4) into (3.3), and take Vh = <fim, with 1 < m < iV — 1, 

N ^ N N 

n—0 n—0 n—0 



(3.5) 



N 

V 

n=0 



^ 1 
^fn(0na;, 4>m.x) = -^{S^ , </>ma;), for 1 < TO < TV - 1. (3.6) 



This set of equations, supplemented with the boundary conditions (3.5) can be 
easily integrated by an ODE solver. 

3.2 The results 

For demonstration purpose, we take as a test case v = 0.2, gi = 0, 52 = and 
h = — sin(^^ + — ). It is easy to check that, for this test case, both the zeroth 
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and first order compatibility conditions at the right corner are met, but those 
at the left corner are not. 

To study the accuracy of the numerical scheme we need a means to mea- 
sure the errors in the solution. Given arbitrary initial and boundary conditions 
in (2.1), generally no analytic solution is known for the Burgers equation, and 
hence there is no way to compute the real errors. As an alternative, we compute 
the comparative errors, which are the differences between two numerical solu- 
tions for the problem, one with the stated mesh sizes, and the other with finer 
mesh sizes. In what follows the term error is to be understood in this sense. 

We first compute the solution of (2.1) without any correction procedure, 
i.e. with S* = in (3.5) and (3.6). The solution is plotted in Fig. 3.2 (a), and 
it displays sharp gradient around the left corner of the time-space axes due 
to the incompatibility between the initial and boundary conditions there. In 
order to have an overview of the structure of the errors in the solution we plot 
the pointwise errors in Fig. 3.2 (b). The pike appear near the left corner of the 
time-space axes, as expected, and it dissipates away quickly. For comparison we 
plot, in Fig. 3.3; the solution and the pointwise errors computed with Correction 
Procedure 1, and, in Fig. 3.4, the solution and the pointwise errors computed 
with Correction procedure 2. With Correction procedure 1 the magnitude of 
the errors at the left corner (see Fig. 3.3 (b)) has been reduced by roughly two 
orders; with Correction procedure 2 the magnitude of the errors at the left corner 
(see Fig. 3.4 (b)) are further reduced. 

The evolution of the maximum errors at each time step is plotted in Fig. 3.5 
(a). The evolution of the maximum errors, with Correction Procedure 1 enabled, 
is displayed in Fig. 3.5 (b). The magnitude of the maximum errors has been 
reduced by roughly two orders (compared to Fig. 3.5 (a)). 

When Correction Procedure 2 is enabled, we see further improvement in the 
accuracy, though it is less dramatic. For comparison we plot the result in the 
same figure as that for the result with Correction Procedure 1, and we see that 
the magnitude of the maximum errors is roughly halved. 

Remark 3.1. In practice a correction procedure should be enabled during a short 
period of time at the beginning, and be disabled afterwards, when the solution 
has become smooth enough, to avoid the computational burden associated with 
the singular terms in (3.3). For computations that produced Fig. 3.5, however, 
we run the simulation for a short initial period of time only, and enable the 
correction procedure for the whole process to avoid unnecessary complications 
in programming. 

Remark 3.2. The zone of rapid variation of the solution and of the error should 
not be confused with a boundary layer effect. Indeed, firstly v = 0.2 is too 
large to produce a sharp boundary layer, as the boundary layer size is of order 
= 0.447 which is essentially of order 1. Furthermore when v is small enough 
to produce a sharp boundary layer, it usually appears for all times, at either 
or both ends of the interval, x = 0,1, whereas in our example the zone of rapid 
variation is concentrated near x = 0, for a short time. Regarding the boundary 
layers for Burgers equation, see e.g. [1, 18], and also [5]. We intend to address 
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Without correction Without correction 




(a) (b) 



Figure 3.2: The solution, and pointwise errors without any correction procedure. 

in a future work the issue of the incompatible data as a boundary layer at time 
t = 0; see as already a first result in this direction in [4]. 

3.3 Comparison of convergence rates 

In this subsection we study the decay of the maximum errors, with and without 
the correction procedures applied. The results also demonstrate the eff'ectiveness 
of the correction procedures. 

Fig. 3.6 (b) shows that, with and without the correction procedure applied, 
the maximum errors at a fixed time t = 0.05 decay at approximately the second 
order (the slope of the line), which is the order of accuracy of the finite clement 
scheme. However, the maximum errors of the simulation with Correction Pro- 
cedure 1 applied is at a magnitude about one order smaller than without any 
correction procedure. The maximum errors of the simulation with Correction 
Procedure 2 applied is smaller by another magnitude of about 0.3. 

The most interesting and informative comparison can be made between the 
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With Correction procedure 1 



With Correction procedure 1 




decay rates of the maximum errors at the initial time step (t = At, At varying 
with different configurations). In Fig. 3.6 (a), the curve for the simulation 
without any correction procedure stick to the upper frame of the figure; the 
maximum errors wont come down whatsoever. With Correction procedure 1, 
the maximum errors decay at roughly the first order with respect to Ax, and 
with Correction procedure 2, the results arc slightly better in terms of magnitude 
of the maximum errors. 



4 A nonlinear reaction-diffusion equation 

In this section we extend the previous correction procedures to the nonlinear 
reaction-diffusion equation: 
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With Correction procedure 2 



With Correction procedure 2 




Ut — vu^rc + p{u) ~ < X < l,t > 

u(0,i) = = 52(0,^ > 0, (4.1) 

u{x, 0) = 

where v > Q\s a. parameter representing viscosity; gi(t),g2(t) and h{x) are given 
functions; and p{u) is a polynomial in u of odd order and with a positive leading 
coefficient. The effectiveness of these correction procedures will be demonstrated 
by numerical results. 

We first derive the compatibility conditions between the initial and boundary 
conditions of (4.1). as explained in Section 2.1: 

0%rder: 5i(0) = /i(0), .92(0) = /i(l), (4.2) 

V'order : 5^(0) = lyh" {0) ~ p{h{0)), g^tiO) = J^h'\l) ~ p{h{l)). (4.3) 
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Without correction 



X 10" 



With correction 




-e — Correction procedure 1 
— Correction procedure 2 



0.01 0.02 0.03 0.04 0.05 
(a) t 




0.01 0.02 0.03 0.04 0.05 
(b) t 



Figure 3.5: Evolution of the maximum errors till t ~ 0.05. (a) Without correc- 
tion procedure, (b) With correction procedure 1 and 2. 



For simplicity, in what follows, we only consider the case where incompatibilities 
are present only at the left corner. Thus we let 

«o = 5i(0)-/i(0)^0, (4.4) 
ai = 5u(0) - iyh"{0) + p{h{0)) ^ 0. (4.5) 
(4.6) 

Of course the method we present for deriving the correction procedures would 
also apply to the incompatibilities at the right corner. 

4.1 The correction procedure 

As we did for the nonlinear convection-diffusion equation in Section 2 and 3, 
we shall combine the correction procedures with the weak formulation of (4.1)i, 
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Decay of the maximum errors at initial time steps 



- Without correction 
-X — Correction procedure 1 
-e — Correction procedure 2 




Decay of the maximum errors at final time steps 

I I I — 

Without correction 
Correction procedure 1 " 
Correction procedure 2 




1 1.5 2 2.5 

logarithm of number of segments in x (base 10) 



1 1.5 2 2.5 

logarithm of number of segments in x (base 1 0) 



Figure 3.6: Decay of the maximum errors, (a) at the initial time step , (b) at 
final time step {t = 0.05). 



for the reason stated in Section 2.2. We multiply (4.1)i by u € 2)(0, 1) and 
integrate by parts to obtain the weak formulation of the equation: 

{ut, v) + i^iu^,v^) + {p{u),v) = 0. (4.7) 



We intend to present the correction procedures (corrections of the incom- 
patibilities to various order) in a unified way. To this end we shall employ the 
following notation: 

{0 if no correction procedure is applied, 

ao'5'o if Correction procedure 1 is applied, (4.8) 

ao'5'o + aiSi if Correction procedure 2 is applied. 

Here the singular corner functions Sq and 5*1 are defined as in (2.3) and (2.11) 
respectively. We let N be the number of segments in the interval [0,l],h — 1/N 
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and let Vh be the finite element space spanned by {ipa, ipi, • • • , (pjv) (see Fig.l). 

We look for the solutions of (4.7) in the form 

uc^ S + Vh{x,t), (4.9) 

where Vh{-,t) e Vh for a.e. t. Plugging (4.9) into (4.7), and taking v — Vh £ Vh 
with w/i(0) ~ Vh{^) — 0, wc obtain 

{vht,vh) + v{vhx,vhx) + {p{vh + S),Vh) 0. (4-10) 
Imposing the boundary conditions and initial conditions in (4.1) wc have 

ivh{Q,t)^g^{t)-S{0,t), 
\vh{l,t)=g2{t)-S{l,t), 

and 

Vh{x,Q)^h{x). (4.12) 
We will solve (4.10), (4.11) and (4.12) for Vh, and then we recover u by (4.9). 

Concerning the compatibility conditions for Vh wc make the following ob- 
servations. For the zeroth order correction procedure, S = aoSo, and by (4.11) 
and (2.4), we have 

(vh{0,t)^gi{t)~ao, ^^^^^ 
I Vh{l, t) = g2{t) - aoS'o(l, i). 

The initial and boundary conditions for Vh satisfy the zeroth order compatibility 
condition. Indeed, 

gi(0)-ao = /i(0) (4.14) 
The compatibility conditions at the right corner are not affected. 

For the first order correction procedure S = a^Sf) + aiSi^ and by (4.11), 
(2.4) and (2.12), 

I Vh{Q, t) = gi{t) - ao - ait, 

y Vh{l,t) = g2{t) - ao<So(l,t) - aiS'i(l,i). 

It is easy to see that the initial and boundary condition for Vh satisfy the zeroth 
order compatibility conditions. They also satisfy the first order compatibility 
condition. To see this, we insert (4.9) into (4.1)i and obtain 

Vht - I'Vhxx + p(ao'S'o + "i-S*! + Vh) = 0. (4.16) 
We first calculate u;it(0,0) from(4.15): 
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Without correction 



Without correction 




vht{0,0) ^ git{0) ~ ai (4.17) 

Then wc calculate the same quantity from (4.16), using the initial condition 
(4.12) instead: 

VhtiO,0) = i^h,,{0)-p{h{0)). (4.18) 
The right-hand sides of (4.17) and (4.18) are equal according to (4.5) 



4.2 The numerical results and interpretation 

In the above we presented the correction procedures for reaction-diffusion equa- 
tion with polynomial reaction terms. For the numerical testing of these correc- 
tion procedures we consider the specific case where 

p{u) = u^. (4.19) 
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With Correction procedure 1 



With Correction procedure 1 












Figure 4.2: The solution, and the pointwise errors with Correction procedure 1 
applied 



As in Section 3.2, we take v — 0.2, = 0, 32 = and h ~ sin(^a:: + jtt). 
It is easy to check that, for this test case, both the zerotli and the first order 
compatibility conditions at the right corner are met, but neither of them are 
met at the left corner. 



Generally, given arbitrary initial and boundary conditions, it is impossible 
to find an analytic solution for the nonlinear reaction-diffusion equation. There- 
fore, in general, it is impossible to compute the real errors. As an alternative, 
we compute the comparative errors, which are the differences between the two 
numerical solutions for the problem, one with the mesh under consideration, 
and the other one with a finer mesh. In the following the term error is to be 
understood in this sense. 
We first compute the solution of (4.1) without any correction procedure, 

i.e. with 5 = in (4.15). The solution is plotted in Fig. 4.1 (a), and it dis- 
plays a sharp gradient around the left corner of the time-space axes due to 
the incompatibility between the initial and boundary conditions there. In or- 
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With Correction procedure 2 With Correction procedure 2 




Figure 4.3: The solution, and the pointwise errors with Correction procedure 2 
apphed 

der to have an overview of the structure of the errors in the solution we plot 
the pointwise errors in Fig. 4.1 (b). The pike appears near the left corner of 
the time-space axes, as expected, and it dissipates away quickly. For compar- 
ison we plot, in Fig. 4.2, the solution and the pointwise errors computed with 
Correction Procedure 1, and in Fig. 4.3, the solution and the pointwise errors 
computed with Correction Procedure 2. With Correction procedure 1 the mag- 
nitude of the errors at the left corner (see Fig. 4.2 (b)) has been reduced by 
roughly two orders (compared to Fig. 4.1 (b))! The most pronounced errors 
actually appear later in the simulation due to accumulation of the errors. With 
Correction procedure 2 the errors (Fig. 4.3) are further reduced by roughly one 
half (compared to Fig. 4.2). 

The evolution of the maximum errors at each time step is plotted in Fig. 4.4 
(a). The evolution of the maximum errors, with Correction Procedure 1 enabled, 
is displayed in Fig. 4.4 (b). The magnitude of the maximum errors has been 
reduced by indeed two orders (compared to Fig. 4.4 (a)). 

When Correction Procedure 2 is enabled, we see further improvement in the 
accuracy, though it is less dramatic. For comparison we plot the result in the 
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Without correction X 1 With correction 
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Figure 4.4: The evolution of the maximum errors tih t ~ 0.05: (a) without 
correction; (b) with Correction procedures 1 and 2. 

same figure (Fig. 4.4 (b)) as that for the result with Correction Procedure 1, 
and we see that the magnitude of the maximum errors is roughly halved. 

4.3 Comparison of the convergence rates 

In this subsection we study the decay of the maximum errors as functions of 
the mesh size(number of the segments in a;), with and without the correction 
procedures applied. The results confirm the effectiveness of the correction pro- 
cedures. Fig. 4.5 (b) shows that, with and without the correction procedure 
applied, the maximum errors at a fixed time t — 0.05 decay at approximately 
the second order (the slope of the line), which is the order of accuracy of the 
finite element scheme. However, the maximum errors of the simulation when 
Correction procedure 1 is applied is about one half order smaller in magnitude 
than without any correction procedure. The maximum errors of the simulation 
with Correction procedure 2 applied is smaller by an additional factor of about 
0.25. 

The most interesting and informative comparison can be made between the 




19 



Decay of the maximum errors at the initial time steps Decay of the maximum errors at the final time steps 
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Figure 4.5: Decay of the maximum errors, (a) at initial step.(b)at final step. 

decay rates of the maximum errors at the initial time step {t = At, At varying 
with different configurations). In Fig. 4.5 (a), the curve for the simulation 
without any correction procedure stick to the upper frame of the figure; the 
maximum errors will not come down whatsoever. With Correction procedure 1, 
the maximum errors decay at roughly the first order with respect to Ax, and 
with Correction procedure 2, the results are slightly better in terms of both the 
magnitude of the maximum errors and the decay rate. 



5 Conclusion 

Incompatibilities between the initial and boundary conditions for evolution 
PDEs have an adverse effect on the accuracy of numerical simulations, especially 
near the time-space corners. No matter how fine the grid is, the magnitude of 
the maximum errors persists, with the pike of the errors moving towards the 
corners as the grid gets finer. 

For the same configuration. Correction procedure 1 reduces the magnitude 
of the maximum errors by about two orders; and Correction procedure 2 further 
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reduces the magnitude by roughly one half. 

For a fixed time, the correction procedures have no effect on the convergence 
rate of the solution at that point, but the correction procedures always give 
more accurate results, with Correction procedure 2 being more effective than 
Correction procedure 1. At the initial time step {t = At), without any correction 
procedure the magnitude of the maximum errors persists as mesh size gets small; 
with Correction procedures 1 or 2, the errors diminish at roughly order 1 with 
respect to Ax, with Correction procedure 2 doing slightly better than Correction 
procedure 1. 

The approach described in this article for deriving correction procedures 
does not depend on any particular property of the viscous Burgers equation or 
the reaction-diffusion equation other than its diffusiveness. Hence we believe 
that these procedures can also be applied to other nonlinear diffusive equations 
in space dimension one. As mentioned earlier, in a work in progress [4] wc study 
a totally different method related for higher space dimensions. 
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